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Abstract. We study the Haldane model with nearest-neighbor interactions. 
This model is physically motivated by the associated ultracold atoms 
implementation. We show that the topological phase of the interacting model 
can be characterized by a physically observable winding number. The robustness 
of this number extends well beyond the topological insulator phase towards 
attractive and repulsive interactions that are comparable to the kinetic energy 
scale of the model. We identify and characterize the relevant phases of the model. 
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1. Introduction 

The integer quantum Hall effect initially [1] and the more extended set of topological 
insulators more recently [2], represent a family of many-body systems with exotic 
transport properties. These properties originate from a robust and hidden topology 
that appears in the bands of the system [3]. More precisely, if such two-dimensional 
non-interacting models are embedded in a torus and diagonalized in momentum space, 
they give rise to a collection of bands, ^^(k), and eigenstates, '0^(k), which are 
characterized by the Chern numbers 

= [dk^ {'(pn\dky1pn) “ (^/’n1V’n)] G (1) 

When the Fermi energy lays in the gap between two bands and the total Chern number 
of the occupied bands is nonzero, the resulting state is a topological insulator. Its key 
physical property is that it is non-conducting in the bulk, but it has robust conducting 
states that live on the boundary of the lattice. These edge states are well known for 
their robustness under perturbations: the nonlocal nature of the topological phase that 
supports them, allows edge states to survive under extreme conditions of distortion, 
impurities and other external perturbations. This resilience of topological edge states 
is highly favorable when considering their realization and detection in the laboratory. 

An important open question in topological insulators is what happens to their 
topological properties in the presence of interactions. One special situation is when we 
force the energy bands of the insulator to be extremely narrow, having 5n(k) almost 
independent of k, and introduce a strong repulsion or attraction. A paradigmatic case 
is the fractional quantum Hall effect (FQHE) [4], where states with fractional filling 
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present strong correlations and excitations with anyonic statistics that could be used 
for quantum computing [5]. 

A complementary situation is when there are obvious dominant energy scales and 
the interaction term is comparable to the kinetic term that originates the topological 
phase. This regime, which has been addressed in earlier theoretical works [6l[7], is 
relevant to what could be experimentally achieved in present quantum simulations 
of topological insulator with ultracold atoms. Such simulations have been recently 
realized in the laboratory for the Harper-Hofstadter Hamiltonian and the 

Haldane model m- 

In this work we employ the Haldane model to study the role of interactions in 
the survival or destruction of the topological phases. This model m is a particular 
instance of a topological insulator or anomalous quantum Hall effect. It is realized 
on a honeycomb lattice with complex hopping amplitudes, as shown in figure We 
rely on a variant of the model introduced in Ref. [12], which is amenable to being 
implemented using ultracold atoms in optical lattices. It was also shown in Ref. m 
that it can introduce new regimes, such as a topological semimetal. This particular 
implementation naturally gives rise to nearest-neighbor interactions [14]. This is 
exactly the setting we simulate here. 

Our study relies on an established tool for characterizing non-interacting 
topological phases, the so called winding number. This is a mathematical object 

1/ = T y d2kS(k) • [VS(k) X VS(k)], (2) 

which can be constructed when the unit cell of the lattice is represented as a pseudospin 
and the eigenstates of the model are represented by unique Bloch vectors, 

S(k) := (cr) = tr(cr |V’o(k)) (V’o(k)|). (3) 

For the non-interacting case the winding number is identical to the Chern number 
m- It has been shown that n is an observable that can be directly measured in 
ultracold atoms experiments [12]. This approach has been generalized to topological 
superconductors [16] and to composite problems (i.e. unit cell dimensionality 2^^), 
where it still signals the existence of topological phase [T5] . In this respect, we regard 
the winding number as a global order parameter in its own right that can be measured 
exp er iment ally. 




Figure 1. (a) Honeycomb lattice model depicting the nearest neighbor, t, next- 
nearest-neighbor hoppings, ta^bi cind the interaction U. (b) Lattice unit cell and 
indices, together with the orientation of the phase in the nearest neighbor hopping. 
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The main result of this work is that the topological phase that originates in 
the non-interacting topological insulator regime extends towards both attractive and 
repulsive interactions, and is at all times deteetable through the winding number (§. 
For very repulsive interactions a charge density wave order (CDW) develops, as already 
indicated by exact diagonalizations in finite lattices with a fixed number of particles 
per site 0 m. For very attractive interactions the lattice fills up, as we work at 
zero chemical potential regime. This phenomenology can be observed with a very 
straightforward mean-field theory that works with the vector field § as an order 
parameter. This theory shows both the survival of the winding number as well as 
the phase transitions into the topological order and into the CDW. These mean-field 
predictions are confirmed by Matrix-Product-States (MPS, also known as DMRG) 
ground state computations with up to 50 lattice sites, a size that doubles what is 
currently achieved through exact diagionalizations. Finally, the destruction of the 
topological insulator that takes place for large attractive interactions is modeled using 
an extension of the mean-field theory by Poletti et al. im, which shows the appearance 
of a superconducting phase in such regimes. 

The structure of this paper is as follows. In section we introduce the form of 
the Haldane model used here. Section introduces a mean-field theory that works 
with the spin field as order parameter and thus gives at all times a proper value 
of the winding number. Section introduces Matrix-Product-States, the variational 
ansatz that we use to get numerical estimates of the ground state properties in finite- 
size problems. Since both our mean field and our DMRG fail for very attractive 
interactions, sectionj^introduces a BGS ansatz that we can use to describe that regime, 
generalizing the ideas in m- Section introduces the main results, describing the 
phase transitions that take place in this model and how they are characterized through 
the winding number, the mean field and DMRG, or the BGS theory. Finally, section 
[^discusses further implications of this work and possible outlooks. 

2. The model 

We work with a variant of the Haldane model that introduces fiux as a phase on the 
nearest neighbor hopping of a honeycomb lattice, as explained in Ref. [12] . The model 
is extended to consider also nearest-neighbor interactions that take place between the 
different sublattices of the Haldane model. Overall, we can write the Hamiltonian in 
the form H = Hq ^ Hu^ with the Haldane model describing the hopping of particles 
in the honeycomb lattice. Introducing pairs of indices, s = (i, j), running over the 
lattice, as shown in Fig. we can write 

Ho= - Yi + H.c.) - ^ e(ata. - bl^) (4) 

vG{vo,vi,V2} s s 

vG{vi,V2,V3} S 

and a nearest-neighbor interaction 

E E jS+V ( 5 ) 

vG{vo,vi,V2} s 

Here, G {(0,0), (1,0), (0,1), (1, —1)} represent displacements between lattice cells; 
tv = texp(i0v), and the intralattice hoppings, ta^bi as well as the nearest-neighbor 
interaction U and a possible lattice imbalance, e, that will be made zero in this work. 
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The hopping part of the Hamiltonian, t and ta^b^ implements a Dirac-type 
Hamiltonian with two distinct singularities in the Brillouin zone. The t term gives 
rise to the kinetic terms of the effective model, while the ta^ H terms implement a 
momentum dependent mass. The combination of both terms can be rewritten as a 
pseudospin model 

Ho = - [ d2ke(k)S(k)-[«'t(k)o-«'(k)], (6) 

JBZ 

where ±^(k) are the energies of the two bands at the given quasimomentum k. At 
half filling, that is one atom per unit cell, the ground state is obtained by placing one 
particle in the lowest band, which is a single particle state with an orientation of the 
pseudospin T(k)’t‘ = dictated by S. 

As explained in Ref. [12], this model could be implemented using two separate 
optical lattices that store atoms in two different internal states, one for sublattice A 
and a different state for sublattice B. The nearest neighbor hopping would then be 
implemented by laser assisted tunneling between internal states, so that the hopping tij 
would carry the phase difference between internal states at their respective positions, 
txy oc exp[k • (x + y)]. In our setup we have chosen one fixed direction for the phase, 
k, shown in figure [^, so that the complex hoppings can be written as 

= (7) 

This approach has the characteristic that the overlap between neighboring sites, which 
makes it possible to have nearest-neighbor assisted tunneling, also allows for a strong 
nearest-neighbor interaction (c.f. Ref. M)- 

3. Mean field theory 
3.1. General idea 

Our approach towards doing a mean-field study of the Haldane model builds on our 
knowledge of the non-interacting solution: in this case, the ground state wavefunction 
is exactly parameterized by the pseudospin orientation S(k) of the fermion that 
populates the k quasimomentum in the Brillouin zone. Consequently, we may write 
the ground state wavefunction at half filling as 

\GS) = 0k |S(k)) , (8) 

where the tensor product is taken over the Hilbert spaces of each quasimomentum 
mode. 

In the interacting case we expect the nearest-neighbor repulsion to have a 
moderate effect on the fermion distribution, more or less preserving the Fermi sea 
and the topological order associated to the vector field S. In other words, we will 
estimate variationally the ground state properties of the interacting model using the 
wavefunction (§ and the orientations of the pseudospins as variational parameters. 

This mean-field approach is quite unique in that we do not work with two 
variational parameters describing the unit cell in position space, but rather employs 
a quasi-continuous normalized vector field S(k) defined on a lattice of N equispaced 
points in the Brillouin zone. Overall we have to work with 2^ degrees of freedom 
representing the orientation of the pseudospins in such sampling. It is important to 
remark that this method becomes exact in the U = 0 limit. 
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Moreover, the mean-field variational wavefunction (§ continuously interpolates 
between the topologically ordered phase, and other phases that are to be found, such 
as the charge density wave (CDW) that we expect in the limit U +oo, where 
atoms localize in either of the sublattices. This variational ansatz, however, does 
not capture other orders, such as a superconducting phase that should arise for large 
negative U ~ — 2|t|, which have to be analyzed using a different ansatz. 


3.2. Spin model 


We rewrite the Haldane model in momentum space, including the interaction. For 
that we introduce the Fourier transformed operators 



s 


( 9 ) 


and similarly for the h operators. Here k is a set of N discrete momentum operators 
that span the first Brillouin zone (BZ), and r is the lattice position. 

The interaction term is written in position space as 


( 10 ) 

Vm S 

with three displacements vo,i ,2 that connect one site to its neighboring cells. The 
Fourier transform of this Hamiltonian becomes 

^ 6j,^6k.43ak,5(ki-k2 + k3-k4)/(ki-k2). (11) 

kl,2,3,4 


The central exponential, summed over r has lead to a (5(-) that enforces the 
conservation of momentum. Finally, we have the weight function 

= ( 12 ) 

Vm 

At this point we notice that Hu can be decomposed into terms that connect two 
momenta and terms that connect four different momenta. Since we are going to use 
a variational wavefunction of the form (|^, the latter terms do not contribute. Our 
Hamiltonian thus reads 

H = - y]e(k)S(k)-<7+^ Y [^p^paVq/(0) + &p^qaVp/(P - q)] -(13) 

k 

Within our variational subspace of one particle per momentum orbital, we have the 
relations 

= ^ (! + <)> flq^q = “^q)’ ^q«q = ' (1^) 

This allows us to rewrite, up to global constants, 

^ = -5l£(k)S(k)-o--^ 5Z<<^q‘/(p-q)-(15) 

k \ p / p^q 

The second term induces charge density wave order, forcing all atoms to be on one 
sublattice or the other, while the third term scatters particles to different momenta. 
Both terms counteract the effect of the magnetic field e{p)S{p) that tries to enforce 
topological order. 



Winding number order in the Haldane model with interaetions 


6 



Figure 2. Matrix Product State structure running through the 2D honeycomb 
lattice. Dashed lines represent the underlying honeycomb lattice, while solid lines 
represent the MPS bond dimensions. 


4. Matrix Product States ground state 
4 . 1 . General idea 

One alternative to confirm the mean-field predictions would be to do some exact 
diagonalization to compute the ground state wavefunction, evaluating the observables 
that build up the winding number. However, this becomes unfeasible for two combined 
reasons. 

The first one is that we need minimum size of the problem to detect the winding 
number. If we perform a simulation with Li x L 2 unit cells in position space, using 
periodic boundary conditions, that would amount to sampling Li x L 2 points in the 
Brillouin zone. A quick numerical calculation revealed that this sampling is insufficient 
to capture the winding number accurately if L < 3. This means that the smallest 
simulation that could conceivable reproduce the non-interacting result would be a 
4x4 cells lattice with 32 sites. 

The second reason is that given that baseline, performing simulations at half¬ 
filling becomes unfeasible using an ordinary computer and algorithms. Already storing 
the ground-state wave function for the 4x4 lattice demands several gigabytes and we 
were not able to do ground state computations without resorting to state-of-the-art 
diagonalization software and a supercomputer with parallelized Lanczos. Note indeed 
that already the largest clusters that are found in the literature [7] at half filling have 
sizes around 24 sites or 3 x 4 cells, which are insufficient for the winding number 
computation. 

A powerful alternative to exact diagionalizations are tensor-network state 
methods. In particular, we have used Matrix-Product-States to represent the ground 
state and a DMRG-type algorithm to compute the minimum energy state within that 
variational ansatz. An MPS is a variational form of a many-body wavefunction that is 
built by contracting tensors (matrices) in a ID scheme. Using the ordering of lattice 
sites and tensors in figure our ground state is written in the form 

mA])= (16) 

ni ...riAT 

where Um ^ {0,1} is the occupation number of the m-th lattice site. Using a rather 
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straightforward optimization procedure based on local updates of the tensors, it is 
possible to minimize the energy, computing the optimal tensors 


A 


= argmin^ 


(^ 1 ^) 


(17) 


4 . 2 . Winding number from MPS 

After computing the ground state we still have to recover the winding number as 
measured in a time-of-flight experiment [12]. In such experiments, the expansion 
of the atoms in the lattice makes them adopt a Fourier transform of their original 
wavefunctions m, 

4 ^ J d3p^(p)e*P'-^tA„(p)t, (18) 

b] ^ j d^pw{p)e^P’^^Mp4, (19) 

where the quasimomentum is mapped to the position of the expanding atoms, 
p = mr/t. In this formula represents the center of the i-th lattice cell in the given 
sublattice, s G {a, b}. The weight w is the Fourier transform of the wavefunction of an 
atom trapped in one lattice site, and for deep enough lattices it can be taken constant 
over a large domain. Finally, N is an overall normalization. 

The previous expression implies that if we measure the spin texture of the 
expanding atoms 

v(p) = (4't(p)(T5'(p)), (20) 

with this texture will be related to the spin texture of the original 

bands in the lattice 

v(p) OC \w{p)\^ . (21) 

m,n 

In particular, when p = k coincides with a quasimomentum in the Brillouin zone, 
v(k) (X S(k), up to factors that drop out when we normalize v to recover S. 

Consistently with the previous reasoning, we have worked with the MPS 
wavefunction computing the vector field 

(22) 

m,n 

where k G tt/M x [—M/2, M/2)®^, is a finite sampling of the Brillouin zone. Using 
this sampling, we then compute the winding number using an accurate formula for 
the solid angle spanned by every three neighboring pseudospins, S = v/|v|, on the 
momentum space lattice [19], 

i^= Yi ^^(Sk,Sp,Sq), (23) 

with the solid angle approximation 

__ (24) 

|a||b||c| + (a • b)|c| + (a • c)|b| + (b • c)|a| ’ 


tan 


U(a, b, c) 


that is valid for small O. 







Winding number order in the Haldane model with interaetions 
4 . 3 . MPS teehnieal remarks 


There are several important technical remarks regarding the MPS simulations. The 
first one regards the use of conserved quantities to restrict the simulation, imposing 
for instance, a total number of particles in the lattice. We have chosen not to do this, 
looking for the ground state of the full Hamiltonian, which is equivalent to minimizing 
the free energy at zero chemical potential. This means that in these simulations 
an increase or decrease of the filling is possible, though, as we will see below, it is 
irrelevant for the topological order as detected by the winding number -a signature of 
the robustness of this quantity, as shown already in [12]. 

The second remark regards the size and type of lattice. MPS simulations were 
done for lattices with 4 x 4, 4 x 5 and 5x5 cells (32, 40 and 50 sites), the latter showing 
the clearest signal. We explicitly use open boundary for the MPS wavefunction because 
of several reasons. First of all, open boundaries are the ones that best reflect a potential 
experiment with ultracold atoms. Second, we wish to show that the winding number 
measurement is a robust and powerful observable that does not restrict to abstract 
geometries, such as tori and spheres. Third, since we already have periodic boundaries 
in the mean field theory, agreement between this theory and the OBC MPS represents 
a strong signature that our approach is correct. 

Finally, we have to comment on the size of the MPS, the so called bond dimension. 
This dimension is the size of the matrices in (16) and it relates to the maximum amount 
of entanglement that is available in a bipartition of the state. The simulations that 
we show here are done with bond dimension y = 100. This restriction is based 
on the need to scan in detail the whole parameter space and the use of long-range 
interactions induced by the 2D-to-lD mapping. However, in this particular study, 
for the observables that we computed, including density correlations and the full 
winding number, we have verified that convergence starts already at very low bond 
dimensions, such as 30 for the 5x5 lattice. We believe that this is possible because 
PEPS can already represent topological insulator phases and topologically ordered 
states faithfully [20] with moderate bond dimensions. Since PEPS may be rewritten 
as MPS with a linear increase in the bond dimension (y ^ y x L, where L is the lattice 
diameter), our MPS representation, which is easier to operate numerically, does not 
need to be too complex to reproduce many of the relevant physical properties. 


5. BCS Ansatz 


We can further analyze our system, in the case of attractive {U < 0) interactions, by 
means of a using a BCS ansatz. In order to do so, and to facilitate of comparison, we 
will follow the procedure used in Ref. m, which deals with a similar tight-binding 
Hamiltonian. Up to trivial factors, the Hamiltonian presented there corresponds to 
the 4> = 0 case in Eq. 0. The BCS ansatz is based on a Wick expansion of the 
interaction term (overwrite term here) in momentum space which gives 

Hint = E (^k^-kflfc + AkOkfoLk - Ak(ak6Lk)) , (25) 

k 


where Ak = -1/A^^ /(^ - k')(&-k'ak') is computed through a self-consistent 

calculation. The order parameters which characterize the phase are nearest-neighbor 
two particle pairing correlators 


S = 


ds^s+vo) : (^s^s+vi) : {dsb> 


'S+V2 


>) 


(26) 



Winding number order in the Haldane model with interaetions 


9 



Figure 3. (a-b) Winding number as a function of the imparted phase and 
the nearest neighbor interaction, t/, using t = l,ta = —tb = 0.1. We plot the 
outcome of the mean-field calculation (a) and of a MPS simulation (b) with 5x5 
unit cells (50 sites), either in momentum (a) or in position space (b). (c) Mean 
field average value of |{cr^}|. (d) MPS expected value of the phase separation, 
{^a,x?T'6,x{ (solid), and of the average number of particles per site, (ria^b) (dash- 
dot), as a function of the interaction, U/t, for various fiuxes ^/tt = 0 (line), 0.5 
(star) and 1.5 (circle). 


where Vm are again the directions spanned by the three nearest neighbors. Following 
Ref. [17], we rewrite the order parameter as ^/|^| = where the are a 

basis for the Ss cyclical permutation group: ui = (1,1, l)/\/3, U 2 = (2,-1, 
and U 3 = (0,1, —l)/v^. A non-zero value of any Wm coefficient (that is, of |^|) 
signals a superfluid phase, while the particular coefficient renders information about 
the symmetry of that phase. As we will see below, the results that we obtain in the 
attractive regime are consistent with those of Poletti et al m- 


6. Results 

Figure [^summarizes the main results from the previous numerical methods. In figures 
[^-b we plot the winding number of the mean-field or MPS simulation of the ground 
state. Both figures show a strong qualitative agreement, exhibiting a phase transition 
from a trivial phase around zero flux, $ = 0, into a topological phase at larger fluxes 
where the Haldane phase is obtained for a wide range of interactions. 

Along the vertical axis we find that the topological phase described by the winding 
number disappears for repulsive interactions around U ^ t. As shown in figure]^, this 
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Figure 4. Results from the BCS ansatz. In figures (b-d) we plot the abs olut e 
value of the order parameters, for k = 1,2,3, from the BCS ansatz ( |26| >, 
while figure (a) shows the the norm |<5|. The solid line delimits where the 
superconducting order parameter appears. Note that this happens well below 
the line U = —1 where our winding number simulations fail. 


transition is due to a spontaneous symmetry breaking of the expectation values, 
signaling the transition into a charge density wave where particles are polarized into 
one sublattice or the other. 

To confirm the phase separation or CDW order, we have studied a similar order 
parameter in the MPS simulation. In particular, we have computed the expectation 
value {uas'^bs)^ which measures the coexistence of atoms in neighboring sites in the 
honeycomb lattice. For large U this value decreases continuously down to zero, 
confirming the separation of species in the lattice. Interestingly, from the point of 
view of the CDW order, this looks like a cross-over, but from the point of view of 
the winding number this looks like sharp phase transition into a disordered regime. 
Whether this is an artifact of poor resolution on the numerical side (for instance 
because the norm of v becomes so small that our computation of the winding number 
is inaccurate), is something that our MPS simulations cannot resolve. 

We have also studied what happens for negative or attractive interactions. When 
< 0, the ground state configuration at zero chemical potential has a filling fraction 
larger than 1/2, that is more than one particle per unit cell. Despite this enlarged 
filling, the ordered phase still persists until U = —1, as evidenced by the winding 
number [Figure |^-b]. At this critical interaction the lattice becomes perfectly 
filled [Figure 1^], forming a Mott insulator with two particles per site. This trivial 
configuration cannot be reproduced with the mean-field, because that wavefunction 
assumed half-filling. In order to study this region of strong attractions we have to use 
the BCS ansatz developed above Figures [^-d show the outcome of that ansatz. 
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revealing that the superconductor phase that was predicted for the honeycomb lattice 
m exists also in the full Haldane model. Most important, the symmetry of that phase 
is strongly dependent on the breaking of time-reversal symmetry in the non-interacting 
limit, i.e. of the topological phase of the U = 0 Hamiltonian. 

7. Discussion 

Summing up, our study reveals that the topological insulator phase survives for a wide 
range of interaction and that this phase is faithfully detected by the winding number 
operator. Given that earlier exact diagionalizations computed the Chern number 
in selected regimes of parameters 1313 , this leads us to conclude that the winding 
number can also reproduce the Chern number in interacting systems with moderate 
correlations. 

The ideas put forward in this work have very straightforward extensions to other 
models, including other topological insulators and composite systems m, and other 
types of interactions, such as on-site Hubbard terms or long-range interactions. In all 
these systems it would be interesting to see whether the winding number may act as 
a precursor of other strongly correlated phases. 

One may also wonder about the applicability of our results to the recent 
beautiful experiment demonstrating the Haldane model in an optical lattice [TO] . This 
experiment uses laser assisted tunneling to implement relying on the original 
lattice to supply t. Compared with our original proposal [12], it has the problem that 
the small overlap between neighboring sites would lead to a small value of U. In this 
case it would be more advantageous to implement other types of interactions, such 
as relying on two-level atoms and implementing on-site repulsion or attraction. Such 
models fall out of the scope of this work, but could be studied using a generalization of 
our MPS and mean-field methods above, combined with the study of partial winding 
numbers m and how they relate to the global topological order. 

Finally, we would like to emphasize the mean-field ansatz which could be of 
interest to other contexts and models. Such momentum-representation mean-field 
theory differs from other mean-field theories that have been developed in position 
space [21] , and connect to long-range interaction classical spin models that have been 
long analyzed in the literature. The connection between these models, the symmetries 
of the underlying topological insulator and how these are broken by the interaction 
terms could be a profitable avenue to understand the nature of the phase transitions 
that have been found in this study. 
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